Morphogen Profiles Can Be Optimised to Buffer Against Noise 
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Morphogen profiles play a vital role in biology by specifying position in embryonic development. 
However, the factors that influence the shape of a morphogen profile remain poorly understood. 
Since morphogens should provide precise positional information, one significant factor is the robust- 
ness of the profile to noise. We compare three classes of morphogen profiles (linear, exponential, 
algebraic) to see which is most precise when subject to both external embryo-to-embryo fluctuations 
and internal fluctuations due to intrinsically random processes such as diffusion. We find that both 
the kinetic parameters and the overall gradient shape (e.g. exponential versus algebraic) can be 
optimised to generate maximally precise positional information. 



PACS numbers: 87.18.-Tt, 87.17.-Pq, 87.10.-e 

Morphogens are signaling molecules which play a vital 
role in biological development by inducing responses in 
a concentration-dependent manner In the standard 
model of morphogen gradients, morphogen proteins orig- 
inate from a localized source, diffuse and are degraded, 
setting up a concentration gradient across the system. 
This gradient can control patterns of gene expression, 
where, for example, a gene is switched on when the con- 
centration is above a certain fixed threshold, but is off 
otherwise. In this work, we focus on morphogen gra- 
dients specifying boundaries of gene expression at fixed 
absolute distances from the morphogen source, a scenario 
that is frequently realised in developmental biology. 

Developmental systems need to be robust to sources 
of noise in order to generate precise patterns of gene ex- 
pression. Here, we address a simple question: in the pres- 
ence of noise, which morphogen profile is most precise in 
specifying the positions of gene expression boundaries? 
In principle, any spatially non-uniform profile could be 
used to position gene expression boundaries; our goal is 
to understand which profiles might be preferred. Previ- 
ous theoretical approaches have predominantly examined 
robustness of morphogen profiles to embryo-to-embryo 
fluctuations in the morphogen production rate 0, S 0] , 
and analysis suggests that algebraic profiles are most pre- 
cise However, even in (hypothetical) embryos with no 
embryo-to-embryo fluctuations, there would still be vari- 
ation in positional information due to internal fluctua- 
tions in morphogen production, diffusion and degrada- 
tion . Such fluctuations impose limits on the precision 
of biochemical signaling @, 0, [1], and could in princi- 
ple alter the shape of the profile with the best precision. 
Internal fiuctuations arc particularly large if the mor- 
phogen is, directly or indirectly, a transcription factor. 
In that case, the arrival of a morphogen molecule at the 
nanometer-scale DNA binding sites on the target genes 
will be a rare, stochastic event [1, Q. Moreover, internal 
fluctuations in the mor pho gen production rate may also 
play an important role [10|. 

Recent experiments in the fruit fly Drosophila 
melanogaster have quantitatively studied the morphogen 



proteins Bicoid (Bed) 0, [U, [13] , Decapentaplegic (Dpp) 
13] and Wingless [ij]. Interestingly, the observed pro- 
files are exponentially decaying in all cases. In order to 
better understand this finding, we compare three classes 
of morphogen profiles (linear, exponential, algebraic) to 
see which is most precise when subjected to the combined 
effects of both external and internal fiuctuations. We 
find that the kinetic parameters describing morphogen 
profiles can be optimised to buffer against the combined 
effects of both sources of noise. By comparing optimised 
profiles, we then see that the overall shape of the pro- 
file (e.g. exponential versus algebraic) can also be opti- 
mised. Exponential profiles frequently emerge as the best 
compromise: such profiles are not particularly robust to 
either external or internal fluctuations taken singly, but 
when both types of fluctuation are taken together, ex- 
ponential profiles can be the most precise. We therefore 
propose a simple design principle for morphogen profiles, 
namely that evolution has selected gradients with op- 
timal robustness to the combined effects of embryo-to- 
embryo and internal noise. Depending on which source 
of noise is most important, qualitatively different mor- 
phogen profiles will be selected. Given that very high 
positional precision can be achieved by morphogens (e.g. 
a few percent of embryo length in the Bed system), it 
seems plausible that optimisation may well be exploited 
by evolution. However, there may still be other con- 
straints on morphogen systems, for example, on the in- 
formation capacity of the signaling network [3, [3] . 

Models. We consider three representative models of 
morphogen gradients: a freely diffusing morphogen with 
a source and a sink [l^; diffusion with linear decay; and 
diffusion with quadratic decay [3] ■ Why study these three 
particular models? Diffusion with linear decay leads to 
exponential morphogen profiles that have been observed 
experimentally. Non-linear decay mechanisms leading to 
algebraic profiles are robust to external morphogen pro- 
duction fluctuations 0] and we study quadratic decay 
(decay via dimerisation) as a representative example [2] ■ 
Bone Morphogenetic Protein (BMP) in Xenopus embryos 
is a possible candidate for a morphogen shaped by such 
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effective non- linear degradation [l7|- Finally, the source- 
sink model generating a linear profile provides insight 
into gradients generated by diffusion with localised degra- 
dation [l6j . This is realisable when the morphogen degra- 
dation factors are tightly localised, such as for retinoic 
acid along the anterio-posterior axis of mice embryos [11] . 

A three-dimensional system is considered with a planar 
source of morphogen at x = that produces a flux of pro- 
teins (J) which diffuse (D) through the system (length 
L). In the absence of fluctuations the morphogen con- 
centration is given by the reaction-diffusion equation 

5tp(x, t) = i^VV(x, t) - ^p(x, tr , (1) 

where ip = n[s^^] {n = 1) and ip = a[^m^s^^] [n = 
2). The boundary conditions are Ddxp{x)\x=vi + J = 
and, for n = 1,2, Ddxp{x)\x=L = 0. The source-sink 
model corresponds to "0 = and a sink at the boundary 
X = L (i.e. p{L) = 0). The steady-state solutions to ^ 
depend only on the linear distance from the source in the 
x-direction: 



Ps-s{x) = {J/D){L-x), 
PHn{x) = {JX/D)e-^^\ 

Pquad{x) = ^(.T+Xo)"^, 



(2) 
(3) 
(4) 



where Ps-s,im.quad correspond to the source-sink, linear 
decay and quadratic decay models respectively and A — 
^/D/JI, a = 6D/a, .TO = {UD'^/Jaf/^ with jSH) valid 
for x,xq,\ <^ L. 

External fluctuations. We first investigate robustness 
to external fluctuations solely in J as this is believed 
to be an experimentally relevant scenario The posi- 
tion Xt where the concentration passes through a thresh- 
old {pt) varies as a result of embryo-to-embryo fluctu- 
ations {6J) in the morphogen production rate. Keep- 
ing the threshold concentration fixed, and expanding 
around xt to leading order, the width due to external 
fluctuations is given hy W ~ 1^p[xt) /\{p' [xt))\- Here, 
Ap{xt) = [{p{xt)'^) — {p{xT)y]^^'^ ((•■■) denotes ensemble 
averaging) are the variations due to external embryo-to- 
embryo fluctuations, and p'{xt) = dxp\x=XT- The widths 
for the three models are Ws-s ~ {L — xt){SJ/ J), Wun ~ 
\{5J/J) and Wq^ad ~ (l/3)a;o(5J/J), to leading order 
in 5J / J . Typically, xt ^ L/2 (as in the Bed controlled 
hunchback (hb) gene expression boundary in Drosophila) 
and Xo < 3A, leading to Wguad < Wun < W^^s- 

Internal fluctuations. Within an embryo internal fluc- 
tuations are also an important source of noise. We con- 
sider morphogen production, diffusion and (if appropri- 
ate) degradation to be stochastic processes. The mor- 
phogen concentration is sampled in a region of linear 
size Aa; corresponding to the size of the binding site 
at target genes. Note that incorporating details of the 
binding process is unlikely to alter our results In- 
ternal fluctuations in particle density within the sam- 
pling volume again cause the position where the gradi- 
ent passes through pT to vary, leading to imprecision 
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FIG. 1: Numerical Results. A: Pquad against x, averaged for 
r = 5000s. Line is fit to ([4]). B: Mean and standard deviation 
of /(x), (error bars from 10 simulations). C: Comparing ujguad 
from simulations with ((T} for J — !{•), J — 2(o) and J = 
5(A) molecules /im^^s^^, with various r and xt. D: coquad 
against non-linear decay probability, averaged for r = 1000s. 
Dashed line is prediction from (O, where a — Pj/{1 + cPj) 
[13 (7 = (Sx)^/St, c = 12.5/im~'^s). Parameters unless 
stated otherwise: J = 1 molecules /im~'^s~^; xt ~ ipm; D = 
Do — Q.&lpm? s~^\ P = 10~^; L — 50/im; lattice spacing 
5x = O.Ol/im; time step 5t — 2.5 x 10~''s. 



in the positional information provided by the gradient. 
The width due to internal fluctuations w is again given 
by a; w IS.p{xt) /\{p' {xt))\ 0, but where Ap is now 
due to internal fluctuations. In the source-sink and lin- 
ear decay models the statistics of the particle number 
n{x) due to internal fluctuations arc Poissonian, since 
both models are linear and morphogen production, dif- 
fusion and de gra dation (if present) are all Poisson pro- 
cesses [1, [13, 120I . We present numerical results below 
demonstrating that the particle number statistics in the 
quadratic decay model are also effectively described by 
a Poisson distribution. This result is not obvious, as 
non-linear decay processes could, in principle, generate 
non-Poisson fluctuations. In all cases, we assume that 
diffusion of morphogen proteins is a purely three dimen- 
sional process (there could be one dimensional sliding 
along DNA but this leads to similar fluctuations as three 
dimensional diffusion [2H). 

Simulations. Stochastic simulations were performed 
on a three-dimensional lattice containing discrete parti- 
cles with particle injection on the plane a; = 0, diffusion 
and (if appropriate) degradation. For the quadratic de- 
cay model, simulations were performed using a range of 
parameter values for J, D, P (probability of two par- 
ticles degrading within a single time step given that 
they occupy the same lattice site) and system size L. In 
Fig. [TJ\_ we demonstrate that concentrations measured 
in our simulations agree well with ^ (the finite size cf- 
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fects at large x do not alter our conclusions [24]). To 
confirm that particle number fluctuations are described 
by Poisson statistics for the quadratic decay model we 
calculated f{x) = [{n{x)^) - {n{x))^ - {n{x))]/ {n{x)) , 
where each simulation was averaged for 5000s (Fig.[T}3). 
In three dimensions diffusion is efficient enough to pre- 
vent the build-up of non-Poissonian correlations result- 
ing from the quadratic decay reactions. This result is 
robust to parameter variations in our simulations (data 
not shown). 

Time/spatial averaging. We now consider the effects of 
time and spatial averaging @, Q, which act to reduce to. 
Time-averaging is performed by the down-stream signal- 
processing network, where the timescale is typically given 
by the lifetimes of the transcripts/proteins of the target 
gene. Over an averaging period t, there can, at most, 
be Nr = r/Tind independent readings of the morphogen 
concentration which reduce the measurement width by a 
factor ~ Nt Intuitively, Tind ~ {^xY / Dq, the typi- 
cal timescale for the sampling volume to empty and then 
be refilled by new protein particles via diffusion Q 
(Do is 'local' diffusion responsible for movement into the 
sampling volume, which may not equal D, the bulk dif- 
fusion). A constant associated with time averaging is 
found numerically (see below). We also include spatial 
averaging, motivated by recent experiments on the Bcd- 
hh signaling pathway [9| . By averaging over a number of 
different nuclei/cells, Ngpat, the effects of internal fluctu- 



ations can be further reduced by a factor ^ J^spat i^P 
to a constant of 0(1), which wc have verified numerically 
does not alter our conclusions). From N^pat ~ CDqt 
where C is a constant which depends on the particular 
arrangement of nuclei/cells. 

Width due to internal fluctuations. For the three mod- 
els we find that the widths lu due to internal fluctuations 
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In Fig. [Tp, we compare ([7]) with simulation results for 
a range of parameter values. We find good agreement 
between the two approaches, where ^2 = 0.56 ± 0.06 is a 
fitting parameter. A similar procedure for the other two 
models yields fco = 0.53 ± 0.07 and fci = 0.60 ± 0.05. 

Optimizing kinetic parameters. Importantly, the un- 
derlying kinetic parameters in the linear and quadratic 
decay models can be optimised to minimise ujun and 
^quad at the threshold position xt- The existence of 
optimal decay rates is a general feature of morphogen 
gradients and occurs for all decay exponents n > |22l |. 
Maximum precision is achieved when A = Q and 
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FIG. 2: Optimising morphogen profiles. e;i„ versus A (solid 
line) and equad versus xo (dashed line) with UHnip) and 
u)quad{A). J — 1 .0 molecules /^m^^ s ^ ^ , 5J/J = 0.07, As = 
3 X 10"^^m, D = 15/im^s"S Do = 0.5^im^s"\ xt = 200/im, 
r = 5 minutes, Nspat = 18. The contributions solely from 
external fluctuations, which are linear in A, xo respectively, 
are omitted for clarity. 



xq = Sxt- For the non- linear decay model we verified 
numerically the existence of such an optimal decay rate, 
see Fig. ID. In general, in the experimentally relevant 
range xt ~ L/2, xq ^ A, we find Ws-s < ujii„ < LOquad- 
Although both the morphogen density and slope affect w, 
the value of the slope is the more important: the source- 
sink model has the steepest, and the quadratic model the 
least steep, profile, thereby generating the above order- 
ing. Comparing this inequality with our earlier result 
Wquad < Wlin < Ws-s, the Ordering is reversed, so that, 
for example, the quadratic decay model performs best on 
external fluctuations but worst on internal fiuctuations. 

The robustness of our three models to the combined 
effects of internal and external fluctuations can now be 
compared. Since internal and external fluctuations are 
statistically independent, the total width is given by 
e 



Vw^ + M^. Wc find that the minimum in eim with 
respect to A (and in Squad with respect to xq) becomes 
more pronounced than is the case with internal noise 
alone (see Fig. 2). As a result, the penalty in terms of 
reduced precision when using non-optimized parameters 
is now more severe than in the internal noise-only case. 
Note that the optimal values of A and xq are reduced 
compared to their values when considering internal noise 
alone, see Fig. [2] 

Optimizing morphogen profile shape. We investigate 
which model gives the most precise positional informa- 
tion when subjected to embryo-to-embryo fluctuations 
(parameterised by SJ/ J) and internal fluctuations (pa- 
rameterised by the averaging time, r). For each SJ/J 
and r, we compute the optimal decay rates for the lin- 
ear and quadratic decay models as above. We then build 
up an effective phase diagram in the r — SJ/J plane, 
determining the most precise morphogen profile for par- 
ticular levels of external and internal noise, sec Fig. [3J\. 
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FIG. 3; A: Effective pfiase diagram for relative position of tfie 
three models considered with varying internal and external 
fluctuations. B: Relative precision of models with varying xt 
in presence of internal fluctuations, with fixed 5J/ J = 0.07. 
Colourmap defined in text, L = 500 fim and other parameters 
as in Fig. [2] Solid lines correspond to eii„ = Es-s (black) and 

Qin = ^quad (bluc) . 



The parameters values used are similar to those found 
from experiments on hb expression in the Drosophila em- 
bryo 1, El, E in. J is kept constant between the 
models, representing fixed resource expenditure. The 
colourmap is defined as (e/i„ — ea)/^iin x 100% where Ca 
is the smaller of Cquad and tg-s- The percentage change 
is negative when the linear decay profile is most precise. 
For small times or small 5J/J, profiles from the sink- 
source model are preferred as this model is best able 
to buffer the dominant internal fiuctuations. For large 
5J / J or large r quadratic decay profiles arc selected, as 
now the quadratic model is now able to best buffer the 
dominant external fiuctuations. At intermediate values 
of 5J I J and t, when W and oj are similar, exponential 
profiles are preferred as they provide the best compro- 
mise between robustness to both internal and external 
noise. From such diagrams we can build up a qualitative 
understanding of why particular morphogen profiles are 
selected depending upon the dominant sources of noise. 
We emphasise that our methodology is general and not 
dependent on the specific parameter values used above. 
For example, \i D = Dq then quadratic decay would be 
favoured in a wider region of the phase space, but there 
would still be conditions when linear decay, or the source- 
sink model, would deliver higher precision. If we choose 
a general decay exponent n, then, for a fixed SJ/ J, we 
find a qualitatively similar picture with small n favoured 
at short averaging times, with larger n selected at longer 
times [2^ . 

Multiple targets. In many developmental systems mor- 
phogens are directly involved in the expression of sev- 
eral target genes at different positions throughout the 
system (e.g. orthodenticle [xt ~ 125^m) [25| and hh 
{xT ~ 240/xm) are both targets of Bed). We consider op- 
timising the morphogen profiles at one expression bound- 
ary and then compare how precisely the profiles deter- 
mine a second boundary at a different position (keeping 



5 J/ J fixed). Using the same optimised morphogen pro- 
files as above we find that exponential profiles, unlike 
the other gradients, lead to the best precision at both 
short and long distances for short times (r < ISmins), 
Fig. [3J3. This flexibility is a potential explanation for the 
widespread use of exponential profiles in developmental 
systems. 

We have applied our analysis to simple (though still 
experimentally motivated) models of gradient formation, 
yielding linear, exponential and algebraic profiles. In- 
triguingly, the best characterized morphogens Bed, Dpp 
and Wingless in Drosophila all have exponential pro- 
files and their decay lengths are sig nificantly less than 
their respective values of xt 0, llSj . This is consistent 
with our conclusion that decay lengths adapt to buffer 
the combined effects of internal and external fiuctuations 
and that exponential profiles perform well when buffer- 
ing combined internal-external fluctuations for a range 
of XT (since all three morphogens have multiple target 
genes). However, more complex processes may also be 
involved in the formation of these gradients, including, 
for example, transcytos is, p re-steady-state measurement 
or mRNA gradients 0, HallSl- I'^ some circumstances, 
morphogen systems are able to scale with embryo length 
[13, Ell- Nevertheless, once these additional effects 
are better characterized, our concepts can still be readily 
applied. Attaining maximal robustness to the combined 
effects of internal and external noise may be a power- 
ful unifying principle in understanding the fundamental 
design of morphogen systems. 

We thank Enrico Coen and Richard Morris for insight- 
ful comments. MH thanks The Royal Society for funding. 
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